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Abstract 

This article is based on a talk presented at a conference "Quantum Annealing and Other Opti- 
mization Methods" held at Kolkata, India on March 2-5, 2005. It will be published in the proceedings 
"Quantum Annealing and Other Optimization Methods" (Springer, Heidelberg) pp. 39-70. 

In the present article, we review the progress in the last two decades of the work on the Suzuki- Trotter 
decomposition, or the exponential product formula. The simplest Suzuki- Trotter decomposition, or the 
well-known Trotter decomposition [1-4] is given by 

where a; is a parameter and A and B are arbitrary operators with some commutation relation [A, B] ^ 0. 
Here the product of the exponential operators on the right-hand side is regarded as an approximate 
decomposition of the exponential operator on the left-hand side with correction terms of the second order 
of X. Mathematicians put Eq. in the form 

^.xA^xB ^ ^x{A^B)\0{x'^) 

and ask what correction terms appear in the exponent of the right-hand side owing to the product in the 
left-hand side. They hence refer to it as an exponential product formula. (The readers should convince 
themselves by the Taylor expansion that the second-order correction in Eq. |^ is the same as that in 
Eq. 121. The higher-order corrections take different forms.) 

We here ask how we can generalize the Trotter formula J^l to decompositions with higher-order 
correction terms. We concentrate on the form 

^x{A+B) _ ^p^xA^'ixB ^^xA^pixB , , , ^pmxB _|_ 0(2;™+!) (3) 

or equivalently 

QPixA^P2xB^P3xA^PixB , ^ ^ qPmxB _ ^x{A+B)+o{x'"'^^) j^^-j 

We adjust the set of the parameters {pi,p2, ■ • ■ ,Pm} so that the correction term may be of the order 
of x™"*'^. We refer to the right-hand side of Eq. 10) as an mth-order approximant in the sense that it is 
correct up to the mth order of x. (See Appendix IXI for another type of the exponential product formula.) 

One of the present authors (M.S.) has studied on the higher-order approximant continually [4-37]. 
The present article mostly reviews his work on the subject. We first show the importance of the ex- 
ponential operator in Sect. Q and the effectiveness of the exponential product formula in Sect. 13 We 
demonstrate the effectiveness in examples of the time-evolution operator in quantum dynamics and the 
symplectic integrator in Hamilton dynamics. Section |31 explains a recursive way of constructing higher- 
order approximants, namely the fractal decomposition. We present in Sect. 0] an application of the fractal 
decomposition to the time-ordered exponential. We finally review in Sect. the quantum analysis, an 
efficient way of computing correction terms of general orders algebraically. We can use the quantum anal- 
ysis for the purpose of finding approximants of an arbitrarily high order by solving a set of simultaneous 
equations where the higher-order correction terms are put to zero. We demonstrate the prescription in 
three examples. We mention in Appendix^a type of the exponential product formula different from the 
form (PJ; it contains exponentials of commutation relations. We give in Appendix IbI a short review on 
the world-line quantum Monte Carlo method with the use of the Trotter approximation 



1 Introduction: Why do we need the exponential product for- 
mula? 



First of all, we discuss as to why we have to treat the exponential operator and why we need an approx- 
imant in order to treat the exponential operator. The exponential operator appears in various fields of 
physics as a formal solution of the differential equation of the form 

p{t)=Mfit), (5) 

where / is a function or a vector and M. is an operator or a matrix. Typical examples are the Schrodinger 
equation 

i^^^Pix,t)=n^pix,t) (6) 

(we put h — 1 here and hereafter), the Hamilton equation 



dt V Qit) J \ q{t) 

(see Eq. ()14|l below) and the diffusion equation with a potential 

^P{x,t) = CP{x,t). (8) 
at 

A solution of Eq. (|3J| is given in the form of the Green's function as 

/(0=G(i;0)/(0)=e*^/(0), (9) 

although it is only a formal solution; obtaining the Green's function G{t\ 0) = e*^ is just as difficult as 
solving the equation (|SJ) in any other way. Another important incident of the exponential operator is the 
partition hmction in equilibrium quantum statistical physics: 

Z = Tre-^^, (10) 

where 7i is a quantum Hamiltonian. 

The exponential operator, however, is hard to compute in many interesting cases. The most straight- 
forward way of computing the exponential operator is to diagonalize the operator Ai. In quantum 
many-body problems, however, the basis of the diagonalized representation is often nontrivial, because we 
are typically interested in the Hamiltonian with two terms or more that are mutually non-commutative; 
for example, the Ising model in a transverse field, 

7^ = -^ J.X'^l-r^af, (11) 

and the Hubbard model, 

n = -tY,J2 (^l^j- + 4-^-) + f^I]".T"a- (12) 

In the first example l|ll|l . the quantization axis of the first term is the spin z axis, while that of the 
second term is the spin x axis. The two terms are therefore mutually non-commutative. In the sec- 
ond example H12() . the first term is diagonalizable in the momentum space, whereas the second term 
is diagonalizable in the coordinate space. In both examples, each term is easily diagonalizable. Since 
one quantization axis is different from the other, the diagonalization of the sum of the terms becomes 
suddenly difficult. 

The same situation arises in chaotic Hamilton dynamics. Consider a classical Hamiltonian 

H{p,q)^K{p) + V{q}, (13) 
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where K(p) is the kmetic term and V{q) is the potential term. The Hamilton equation is expressed in 
the form 



di V J \ ^~K{p) )-\k- )\q 
where the operators K- and V ■ are symbolic ones standing for the operations 

k-p=4^K{p) and V-q=4^V(q]. (15) 
dp dq 

Although each operation of K- and V- is simple enough, the "Hamiltonian" operator 

is not easily tractable. This is because the kinetic part and the potential part, 

^ ^ ( ^. ) and V ^ '^-y (17) 

do not commute with each other; see an example in Sect. 12.21 below. 

To summarize this section, we frequently encounter the situation where the exponential operator of 
each term, e^^ and e^^, is easily obtained and yet the desired exponential operator e^^'^^^) is hard to 
come. This is the situation where the Trotter decomposition becomes useful. 



2 Why is the exponential product formula a good approximant? 

We discussed in the previous section the importance of the exponential operator and the necessity of a 
way of treating it. We here discuss a remarkable advantage of the Trotter approximant to the exponential 
operator. 

Let us first confirm that the Trotter approximant is indeed a first-order approximant. By expanding 
the both sides of Eq. (Q), we have 



„xA„xB 

e e 



I + x{A + B) + ^x\A + Bf + 0(a;3) 

I + x{A + B) + {A^ + AB + BA + B^) + 0{x^), 



= \I + xA + -x^A^ + Oix-") \I + xB + -x^B^ + 0(a;-^) 



I + x{A + B) + -x^ {A^ + 2AB + B^) + Oix^ 



(18) 



(19) 



where / is the identity operator. The difference between the two comes from the fact that in the approx- 
imant (|19|l . the operator A always comes on the left of the operator B. Hence we obtain 



qxA^xB ^ ^x{A+B) + ^x'^lA,B]+oix'') 



(20) 



In the actual application of the approximant, we divide the parameter x into n slices in the form 



^(A+B) + i(f)'[A,S]+o((f )') 



x(A+B)A 



-[A,B]+o{ 



(21) 



Thus the correction term vanishes in the limit n — > oo. We refer to the integer n as the Trotter number. 

Now we discuss as to why we should be interested in generalizing the Trotter approximation. The 
Trotter approximant and the generalized one in fact, have a remarkable advantage over other 
approximants such as the frequently used one 



^x{A+B) =i + x{A + B) + Oix''). 



(22) 
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The approximant of the form JJJ) conserves an important symmetry of the system in problems of quantum 
dynamics and Hamilton dynamics. 

In problems of quantum dynamics, the exponential operator, or the Green's function e~'*^ is a 
unitary operator; hence the norm of the wave function does not change, which corresponds to the charge 
conservation. We here emphasize that the exponential product 

g-itpiAg-itp2-Bg-itp3-4 , , , g-itpM-B ^23) 

is also a unitary operator. The perturbational approximant (|22|l . on the other hand, does not conserve 
the norm of the wave function; in fact, the norm typically increases monotonically as the time passes as 
we demonstrate in Sect. 12 .11 below. 

In problems of Hamilton dynamics, the time evolution of the Hamilton system conserves the volume 
in the phase space {p,q\, which is called the symplecticity in mathematics. The exponential product 
formula, in general, also has the symplecticity. 

The time evolution of the Hamilton equation H14|l is described by the exponential operator 

*' f V (24) 



m J V ?tti) 

where Ti is the "Hamiltonian" operator . The Trotter decomposition approximates the time evolution 
with the operator 

e*«~ (e^Kg^v^" (25) 

with K. and V given by Eq. (|17|) . The operator C"'^ describes the time evolution over the time slice t/n 
of a Hamilton system with only the kinetic energy K{p). It thereby conserves the phase-space volume, 
so does the operator e^^. The whole Trotter approximant therefore conserves the phase-space volume. 
This holds for any exponential product formula in the form (O as well. Hence the exponential product 
formula, when used in the Hamilton dynamics, is sometimes called a symplectic integrator. 

In equilibrium quantum statistical physics, the operator e~^'^ does not have a particular symmetry 
except the symmetries of the Hamiltonian itself. The above advantage of the exponential product formula 
is hence lost when applied to numerical calculations of the partition function Z = Tre~^^. In fact, in 
applying the higher-order decomposition (^J to the world-line quantum Monte Carlo simulation, some 
of the parameters {pi,p2, ■ ■ ■ ,Pm} are negative, which causes the negative-sign problem in systems that 
usually do not have the negative-sign problem [38]. The negative-sign problem is the problem that the 
Boltzmann weight of the system to be simulated becomes negative for some configurations. 

Thanks to a recent development of the world- line quantum Monte Carlo simulation [39], the higher- 
order decomposition is not necessary anymore in some cases; the simulation is carried out in the limit 
n — > oo from the very beginning and hence the order of the correction term does not matter in such cases. 
See Appendix for a brief review over the recent development. 

2.1 Example: spin precession 

The fact that the exponential product formula keeps the symmetry of the system is one of its remarkable 
advantages. In the present and next subsections, we demonstrate that this indeed affects numerical 
accuracy strongly. In the present subsection, we use a simple example of quantum dynamics, namely the 
spin precession. 

Consider the simple Hamiltonian 

n^a,+r<j,= ( ^ ^^]. (26) 



If we start the dynamics from the up-spin state 



^(0) - ( J ) ' (27) 

the spin precesses around the axis of the magnetic field H — (r,0, 1) with the period 

T = , 28) 
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Figure 1 : The energy deviation due to the approximations given by (a) the Trotter approximant H29|) and 
(b) the perturbational approximant (|30|l . In both calculations, we put T = 3/4 and At = 0.0001. The 
initial state is the one in Eq. (|27|1 with the energy expectation {H) = 1. 



Although it is easy to compute the dynamics exactly, we here use the Trotter approximant 

G{t + At; i) ~ Q-'^^fy^Q-'^^t-Ca. (29) 

and the perturbational approximant 

G{t + At;t)~I- lAtn = I - iAt{(7^ + Ta^). (30) 

The exact dynamics should conserve the energy expectation (7i). Figure. ^shows the energy deviation due 
to the approximations. The error in the energy of the Trotter approximation (|29|l oscillates periodically 
and never increases beyond the oscillation amplitude. The period of the oscillation in Fig.^a) is equal to 
that of the spin precession. We can understand this as follows: when the spin comes back to the original 
position after one cycle of the precession, it comes back accurately to the initial state H27|l because of the 
unitarity of the Trotter approximation, and hence the oscillation. 

In contrast, the error in the energy monotonically grows in the case of the perturbational approximant 
as is shown in Fig. ^b). This is because the norm of the wave vector increases by the factor 

II l-iZ\m ||~ 1 + Z\t II 7^ ||> 1. (31) 

The remarkable difference between Fig. QJa) and Fig. H^b) thus comes from the fact that the Trotter 
approximant is unitary. 



2.2 Example: symplectic integrator 

We next demonstrate the Trotter decomposition (|25|l in an interesting example of chaotic dynamics. We 
again emphasize that keeping the symplecticity of the Hamilton dynamics has an important effect on 
numerical accuracy. 

Let us first notice that the operators in Eq. (|17|l satisfy 



/C^ = = 0. (32) 



We therefore have 



e^^' = + !0 = L A.^rr..^ ) , (33) 



q \ 1 J \ 1 



e-(^:).(Z + V^.)(f).(^:-^^^-^(^^). (34) 

Note that applying the two operators in the order Q^^^Q^^t is different from applying them in the order 
gVzitgK:zit. jj-^ ^]^g former, the update of g in the application of e'^"^* is done under the updated p, whereas 
in the latter, it is done under p before the update. 
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Figure 2: Simulations of the system 1^^5|l . 
The initial condition is pi = P2 — 0, 
qi =2 and q2 = I with the energy E — 2. 
The time slice is Z\i = 0.0001. (a) The 
movement of the system in the coordi- 
nate space (91,92) for 700 << < 900. 
The broken curves indicate the hyperbo- 
las |(?i(?2| — 2. (b) The energy fluctuation 
due to the Trotter approximation H25fl . 
We plotted a dot every 1,000 steps, (c) 
The energy increase due to the approxi- 
mant (jSEl- 





8 - 
7 - 




6 - 




5 - 






D 


4 - 


c 




w 


3 - 




2 - 




1 - 




- 









200 



Umeno and Suzuki [11, 12] demonstrated the use of symplectic integrators for chaotic dynamics of the 
system 

K{p}^l{pi^+P2^) and V{q)^^qi^q2^. (35) 

The equipotential contour is given by 19192! =constant; hence the system is confined in the area sur- 
rounded by four hyperbolas as exemplified in Fig. EJa) . The exact dynamics should conserve the energy. 
The Trotter approximation of the dynamics, Eq. (|25|l . gives the energy fiuctuation shown in Fig. I^b). 
The energy, though deviates from the correct value sometimes, comes back after the deviation. In fact, 
the deviation occurs when the system goes into one of the four narrow valleys of the potential; it is 
suppressed again and again when the system comes back to the central area. 

This is in striking contrast to the update due to the perturbational approximant 



(/ + Atn) 




^t4-^V{q) 



(36) 



which yields the monotonic energy increase shown in Fig. |2tc) . The reason of the difference between the 
approximants, though less apparent than in the case of the previous subsection, must be keeping the 
symplecticity, or the conservation of the phase-space volume. 



3 Fractal decomposition 

We emphasized in the previous section the importance of the exponential product formula. In the present 
section, we describe a way of constructing higher-order exponential product formulas recursively [5-14] . 
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The easiest improvement of the Trotter formula Q is the symmetrization: 



^2(2;) = ei-^e^^et^ = e-(^+^)+-'^«3+^'^fl5+-, (37) 

The symmetrized approximant has the property 

S2{x)S2i-x) = e5^e"^et^e-t-4e— Sg-f A ^ (3g) 

because of which the even-order terms vanish in the exponent of the right-hand side of Eq. H37|) . We can 
thereby promote the approximant (|37|l to a second-order approximant. 

Now we introduce a way of constructing a symmetrized fourth-order approximant from the sym- 
metrized second-order approximant H37|l . Consider a product 

S{x) EE S2{sx)S2{{l-2s)x)S2{sx) (39) 

_ ^^xA^sxB^^^xA^(l-2s)xB^i^xA^sxB^^xA ^^q-j 

where s is an arbitrary real number for the moment. The expression H37|) is followed by 

S{x) = S2isx)S2iil~2s)x)S2{sx) 

^ ^sx(A+B)+s^x^F/.3+0{x'^)^{l~2s)x(A+B) + {l-2sfx'^R3+0(.x'^)QSx{A+B)+s'-^x^R3+0(x^) 

^ ^x{A+B) + [2s^ + {l-2sf]R3+0{x'') (4J^-) 

(The readers should convince themselves by the Taylor expansion that the third-order correction in the 
exponent of the last line is just the sum of the third-order corrections in the exponents of the second 
line. This is not true for higher-order corrections.) Note that we arranged the parameters in the form 
{s, 1 — 2s, s} in Eq. 139|l so that (i) the first-order term in the exponent of the last line of Eq. (|41|l 
should become x{A + B) and (ii) the whole product S{x) should be symmetrized, or should satisfy 
S{x)S{—x) = I. Because of the second property, the even-order corrections vanish in the exponent of 
the last line of Eq. H41() . Making the parameter s a solution of the equation 

2s^ + 2sf =0, or s = 5—= = 1.351207191959657 • • • , (42) 

2 — V 2 

we promote the product (|39|) to a fourth-order approximant [5] . 

Following the same line of thought, we come up with another fourth-order approximant [5] in the form 



Si{x) = S2{s2xfS2{{l-^S2)x)S2{s2xf (43) 

_ Q^XA^82XB ^82XA^S2XB 'l'^ XA^{1-AS2)XB ^ ^ 2'^ XA^S2XB ^S2XA^S2XB ^^xA (44) 

where the parameter S2 is a solution of the equation 

4s2^ + (1 - 4s2)^ 0, or S2 = = 0.414490771794375 ••• . (45) 

4 — V 4 

We can compare the fourth-order approximants H39|) and H43|) using the following diagram. Suppose 
that the exponential operator is a time-evolution operator from the time t — Q to the time 

t = x. In the product Ip!?^ . the term 52 (sx) on the right approximates the time evolution from t = to 
i = sa; ~ 1.35a;, the term 5'2((1 — 2s)a;) in the middle approximates the time evolution from t — sx to 
t — sx -\- {\ — 2s)x = (1 — s)x ~ —0.35a;, and the term S'2(sa;) on the left approximates the time evolution 
from t = (1 — s)x to t — [1 — s)x + sx = x. Let us express this time evolution as in Fig. OJa). The 
product (|43|l is similarly represented as in Fig. I^Jb). 

As is evident, the first product (|39l) has a part that goes into the "past," or i < 0. This can be 
problematic in some situations; in the diffusion from a delta-peak distribution, for example, there exists 
no "past" of the initial delta peak. The second product (|43|) does not have the problem and hence is 
recommended for general use. 

Once we know how to construct the fourth-order approximant from the second-order approximant, the 
rest is quite straightforward [5] . Following the construction (|43() , we construct the sixth-order approximant 
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Figure 3: Diagrams that rep- 
resent the time evolution of 
(a) the fourth-order approxi- 
mant (|39() and (b) the fourth- 
order approximant H43|l . 



Figure 4: Diagrams that represent the time evolution of (a) the 
six-order approximant H46|l and (b) the eighth-order approxi- 
mant (|^ . 



in the form 



with 



Sq{x) = S4{S4X)'^ S4{{1 — 4:S4)x)S4{s4x)'^ 

= (5'2(S4S22^)^5'2(S4(1 - AS2)X)S2{S4S2X)'^Y 

xS2{{l - As4)s2x)^S2{{l - 4s4)(l - 4s2)a;)52((l - 4s4)s2:e)' 
X (5'2(s4S2a;)^5'2(s4(l - 4s2)x)S2{s4S2x)'^y 

1 



4s5 



(1 - 4s4 



0, 



S4 



and further construct the eighth-order approximant in the form 

Ss{x) EE Se{s6xfS6{{l - Ase)x)Seis6x)'^ 

with 



0.373065827733272 - •• 



(46) 
(47) 

(48) 

4-^ — — • ^''^ 
These approximants are represented by the diagrams in Fig. ^ We can continue this recursive procedure, 
ending up with the exact time evolution, where the diagram ultimately becomes a fractal object. This is 
why the series of the approximants is called the fractal decomposition. It is an interesting thought that 
the back-and-forth time evolution in a fractal way reproduces the exact time evolution. 



44 + (1 - 4se 



0, 



S6 



1 



0.359584649349992 • 



4 Time-ordered exponential 

Before going into another way of constructing higher-order exponential product formulas, let us intro- 
duce, as an interlude, an important application of the exponential product formula. We show how to 
approximate the time-ordered exponential [10]. 

We have considered until now only the case where the operators A and B do not depend on a;, or 
in other words, only the time evolution of a time-independent Hamiltonian. The fractal decomposition 
introduced in the previous section needs modification when applied to problems such as the quantum 
dynamics of a time-dependent Hamiltonian; in quantum annealing [40-42], for example, the transverse 
field r in the Hamiltonian Hll|) is changed in time. 

The time-evolution operator of the quantum Hamiltonian 

n{t) = A{t) + B{t) (50) 
is not simply e^'^* but a time-ordered exponential in the form 



G{t2;h) = T 



exp H{s)ds 



(51) 
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It is quite well-known that 

Gi{t + At-t) = e-^^tA(t+At)^-iAtB{t+At) 

is an approximant of the first order of At and 

G2{t + At-t) = e-^-4*'4(t+izlt)g-izltS(t+i/it)g-^zltA(t+izlt) 



(52) 



(53) 



is an approximant of the second order. How do we construct higher-order approximants? We here show 
that a slight modification of the fractal decomposition gives the answer. 
The key is to introduce a shift-time operator [10] defined in 

F(i)e-'^*^G(t) = F{t + At)G{t). (54) 

Note that the operator acts on the function on the left. The shift-time operator is expressed in the form 



(55) 



in the case where F{t) is an analytic function, but the definition H54|l does not limit its use to the analytic 
case. If we have two shift-time operators, the result is 

F{t)e-'"^''^ G{t)c-'"^''^ H{t) = F{t + At)G{t)e-'"^''^ H{t) 

= F{t + 2At)G{t + At)H{t). (56) 

With the use of the shift-time operator, the time-ordered exponential 151|) is transformed [10] as 

rt+At ^ 



T 



exp — 1 



n{s)ds 



-iAt(H{t)+T) 



(57) 



We can prove this by using the Trotter approximation as follows: 

n — >oo V / 



hm e-^«(*)e-^^e-^«(*)e-^ 



lim e"'^^(*+'^*'e"'^^(*+'^'^*) • • • e"''^^(*+"'^*) 

rt+At ' 



T 



exp — 1 



nis)ds 



(58) 



Decomposing the Hamiltonian into two parts as in Eq. H50|l . we have now three parts in the exponent 
of the time-evolution operator as in 



exp — 1 



t+At 



n{s)ds 



^ ^-iAtiAit)+B{t)+T) ^ 



(59) 



We then approximate the exponential in the right-hand side of Eq. H59|l . The first-order approximant is 
given by 



Gi{t + At]t) 



^-iAtA(t) ^-~iAtB{t) ^-lAtT 
^-iAtA{t+At)^-iAtB(t+At) 



(60) 



and the second-order approximant is given by 
G2{t + At;t) 



^-^Atr^-^AtA{t)^-iAtB{t)^-^AtAit)^-^AtT 



e 2 



^AtA(^t+^At) -iAtB[t+^At) -^AtA(^t+^At) 



(61) 
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Higher-order approximants are given by the fractal decomposition of the three parts, A, B, and T. The 
fractal decomposition of three parts is easily obtained by substituting 



S2{x) = ef '^et^e^'^et^et^ = e-(^+B+c)+o(-^) (62) 
for Eq. (|37|l . The fourth-order approximant is thereby obtained [10] as 

G4{t + At-.t) = (^-h^2AtT^-^S2AtA{t)^~is2AtB{t)^-~-^S2AtA{t)^~^S2AtT\j^ 

y^Q-^{l-4s2) Atr ^~^{l~4s2) At A{t)^-i(l-4s2)AtB(t)^-^{l-is2) At A(t)^-^{l-As2)AtT 
^ .g 2S2AtT^-^S2AtA(t)^-^is2AtB{t)^-^S2AtA{t)^-^S2AtT'^'^ 



^-■^S2AtA[t+?^^At)^-is2AtB[t+?^^At)^~-^S2AtA[t+^^^At) 
^g-^S2AtA(t-|- ^~|°^ At) ^-is2AtB{t+ ^'2^^ At) ^-^S2AtA{t+ ^'2''^ At) 
^^-■^S2AtA[t+^At)^-iS2AtB(t+^At)^-^S2AtA(t+^At) 
^^-^S2AtA{t+^At)^-is2AtB{t+^At)^-^S2AtA[t+'^At) 

^^-^S2AtA(t+^At)^-is2AtB(t+^At)^-^S2AtA[t+^At) (gg-j 



with the coefficient S2 given by Eq. 1)45(1 . 



5 Quantum analysis — Towards the construction of general de- 
compositions — 

In the last section before the summary, we discuss the calculus of the correction terms. In the fractal 
decomposition, we construct higher-order approximants recursively. Is it possible to construct higher- 
order approximants directly, not recursively? In fact, Ruth [43] found (not systematically) a third-order 
formula 

eM^^e^^-^e3^'^e"^^-^e"i^'^e^^ ^ g2:(A+s)-i-o(2;'')^ (64-) 

which would not be found within the framework of the fractal decomposition. 

For the purpose of finding higher-order formulas directly, we need to compute the correction terms in 
the exponent as 

^PixA^P2xB ^P3xA^P4,xB ^ ^ ^ ^pmxB _ ^x{A+B)+x^R2+x^R3 + --- {Q5) 

This is one of the aims of the quantum analysis developed by one of the present authors (M.S.) [29-35]. 
Then we can put the correction terms to zero up to a desired order and solve the set of non-linear 
simultaneous equations 

i?2 = 0, i?3 = 0, • • • , i?™ = 0, (66) 
thereby obtaining the parameters {pi}. 

5.1 Operator differential 

The main feature of the quantum analysis is to introduce operator differential. In order to motivate the 
readers, suppose that we can write down an identity 

ffAf d/(A) dA{x) 
^/(A(x)) = ■ (67) 

where f{A) is an operator functional. The derivative with respect to x on the right-hand side is well- 
defined; for example, dA{x)/dx = B + 2xC for A{x) = xB + x^C. Now, is it possible to define the 
differentiation df{A)/dA? 

Let us discuss as to what should be the definition of the operator differential in order for the iden- 
tity H67f) to hold. The definition of the x derivative is expressed as 

A{x + h)^A{x) + h^^ + 0{h^)- (68) 
da; 
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The left-hand side of the identity (|67|) is given by the definition of the derivative as 

da;'^ ^ ^ h^Q h h^o h ^ ' 

The identity (|67|l suggests that the operator differential d/(A)/dyl must be a hyperoperator that maps 
the operator dA(a;)/da; to the operator given by Eq. 

Thus we arrive at the definition of the operator differential within the framework of the quantum 
analysis [29] : if we can express the operator given by 

d/(A).limZ(A±Z^^^^W(A) 

in terms of a hyperoperator mapping from an arbitrary operator AA as in AA — > d/(j4), then we refer 
to the hyperoperator as an operator differential d/(A)/dA and denote it in the form 

d/(A) = ^-dA (71) 

We stress here that the operator differential d/(yl)/dA must be expressed in terms of A and the commu- 
tation relation of A, or the "inner derivation" 

8a^\A, ], (72) 

but not in terms of the arbitrary operator dA. The convergence of Eq. H70() is in the sense of the norm 
convergence which is uniform with respect to the arbitrary operator dA. 

Let us consider the application of the above in a simple example /(A) = . The definition H7()|l is 
followed by 

,,,,, , {A^hAAf-A^ , HAAA^hAAA^h^iAAf 
df(A) = \va\- = hm — 

h^o h h^O h 

= AdA + dAA^2AdA- {AdA- dAA) 

= {2A - 6a) dA. (73) 



Thus we have [29] 



"^^^'^ -2A- 6a. (74) 



dA 

li A ^ xB + x'^C, we use the result {731) for Eq. (EZl and have 

^ (xB + x^Cf = (2xB + 2x^C- 6^b+x-^c) (B + 2xC) 
dx 

= {2xB + 2x^C) {B + 2xC) - {xB + x^C,B + 2xC] 

= 2xB^ + Ax^BC + 2x^CB + Ax^C"^ - 2x\BC - CB) - x^{CB - BC) 

= 2xB'^ +?,x'^BC + ?,x'^CB ^Ax^C"^, (75) 

which is indeed identical to the result of straightforward algebra. 

We cannot see in this simple example any merit of the use of the quantum analysis. The readers 
should wait for more complicated examples given later in Sec. 15.31 where we show that the differential of 
exponential operators is given in terms of the inner derivation. The Lie algebra is defined by commutation 
relations, or the inner derivation; it is hence essential to obtain results in terms of the inner derivation, 
not in terms of naive expansions such as the right-hand side of Eq. (|75|l . 



5.2 Inner derivation 

We here provide some of the important formulas of the inner derivation (|72|l as preparation for the next 
subsection, where we give the differential of exponential operators. 

First, we have linearity: for any c- numbers a and h, the inner derivation of the operators A and B 
satisfies 

6aA+bB ^[aA + bB, ] = a[A, ]+b[B, ]^a5A + b6B. (76) 
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Any powers of the operator A are commutable with the inner derivation of any powers of the same 
operator: 

[A",<5a"]=0, (77) 

because 

A^Sa'^B = A"" [v4", = [^", A^B] = Sa'^A'^B (78) 

for an arbitrary operator B and any integers m and n. We can generaUze the identity (|77|l to the case of 
any analytic functions of the operator A: 

[f{A),Sg^A)]^0, (79) 
where f{A) and g{A) are defined by the Taylor expansion as 

oo oo 

/(A) = ^a„A" and 5(A) = ^6„^". (80) 

ri=0 n=0 

Next, we prove the identity [29] 

^f(A)a(A) = f{A)6g^A) + 9{A)Sf(^A) - Sg(A)Sf(A)- (81) 
The proof is as follows: for an arbitrary operator B, we have 

fiA)5g^A)B + g{A)Sj^A)B ~ Sg^A)SfiA)B = /(A) [g(A), B] + g{A) [f{A),B] - [g{A), [f{A),B]] 

^ f{A) [g{A),B] + [/(A), B] g{A) = [f{A)g{A),B] 

= ^f{A)g(A)B. (82) 

Note that we can rewrite the identity (|81|l as 

h{A)giA) = Sg(^A)fiA) + g{A)Sf(^A) - Sg(A)Sf(A) 

= Sg(^A)(.fiA)-6f(^A))+g{A)6f(^A) (83) 

because of the identity H79|l . In the special case f{A) — A, we have 

SAg(A) = Sg(^A) {A - Sa) + g{A)5A. (84) 

With the repeated use of the identity 184|l , we then prove the identity [29] 

Sa^ - A" -{A- SaT (85) 

for any integer n. This is proved by means of mathematical induction. The identity (|85|l indeed holds 
for n = 1. Assume now that 

6A,^-^=A''-' -iA-SAr~\ (86) 

Then the identity yields 

Sa^ = SAA'^-^^SA^-i{A~SA) + A"'HA^[A''-'-{A~SAr'^]{A-SA)+A'''^SA 

= A^-iA^SAT- (87) 
An interesting and quite well-known identity is 

^xAq^-xA ^^xSaq^ (88) 

We can prove this by differentiating the left-hand side by x. First, note that 

^^e^^Be^^"^ = e^^ASe-^-^ - e^-^SAe"^-^ = e="^ [A, B] e"^"^. (89) 
da; 

We thereby have the following in each order of a;: 

= e^[A,B]e-^\^^^5AB, (90) 



^e'^Be- 
ax 



= e-^[AM.B]]e-^\^^ = 5A^B, (91) 

x=0 
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which proves the identity (|88|l . As a corollary, we obtain the following identity: 

giAg^B^g^* if e^e-^^e*. (92) 
The proof is straightforward; for an arbitrary operator C, we have 

e^A^Ssc ^ e^e^Ce-^e"^ = e*Ce-* = e**C. (93) 

5.3 Differential of exponential operators 

We are now in a position to discuss the differential of exponential operators. We begin with the differential 
of the power of an operator, f{A) = A^, a generalization of the identity l|74() . The result is [29] 

AA 5a ^ ' 

An important comment is in order. The identity H94|l does not claim that the inverse of 5 a is well-defined. 
In fact, the inner derivation 5a in the denominator is canceled when we expand the numerator of the 
second expression. The denominator is well-defined only in such cases. 

We use the identity H85(l in the derivation of the identity l|94|l . The definition l|70(l is followed by 

df{A) = limIA±i^^^ti!=^A^-i(dA)A"-^ 

nA"-^ - ^ A^-^5Ar.-}j dA = !^nA"^ - ^^^^ ^""^ - - ^aT^^ | dA 

= '^l^ipMliA='-fliA. (95) 

OA OA 

Note again that the transformation in the fourth line is well-defined only because the expansion of the 
numerator cancels the denominator. 

We can generalize the identity (|94|l to any analytic functions defined by the Taylor expansion l|8U|) . 
The result is 

df{A) f{A)-f{A-5A) 5j^ 
dA - 5a ~ 5a- ^ ' 

It is interesting to note that the operator differential or the quantum derivative [10] is expressed by a 
difference form of hyperoperators. As a special case, we arrive at the identity [29] 



e r • (97) 



dA 5a 5a 

5.4 Example: Baker-Campbell-Hausdorff formula 

We now use the formula (|97|l for the derivation of the Baker-Campbell-Hausdorff formula, or the derivation 
of higher-order terms of the exponent '^(x) given in 

e*(")=e"V^. (98) 

The differential of the left-hand side of Eq. (|98|) gives 

dx d<I> dx <5<i.(2;) da; 

owing to Eq. (|97|l , while the differential of the right-hand side of Eq. (|98|l gives 



_d_ 

dx 



^xA^xB ^ ^xA^^xB _^ ^xA^xBg ^ ^xA^xB (g-^B^g^^B + B) = e*'^' (e"^*^ A + B) , (100) 
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where we have used the identity H88|) . Equating the both sides, we have 

d$(x) ^ ^^(.) u-.s.^ fH^) ^ ^...5) (101) 

dx 1 - e~**<"^) c'^'>'(==) - 1 ^ ^ 

The second equahty is due to the identity (|92|l . 

We can expand the right-hand side of Eq. (|101|) as foUows. Note here that 

e-s^cx) = ^xSa^xSb yjgids ,5<(,(^) = log (e^-^-^e^"^^) . (102) 
Thus we transform Eq. pUl|l as 

d£(x^ ^ fog (e-^^e-^^) ^ ^ (-1)^ , _ 1)'= + e^'-B) . (103) 

fe=0 

We finally arrive [30] at 

\k px 



fc=0 ^" -'^ "'0 



(104) 



It is very important to notice here that all the expansion terms are given by commutation relations. One 
of the merits of the quantum analysis is to be able to express the expansion in terms of commutation 
relations. 

Let us derive, for example, the term of the third order of or the second order of t of Eq. (|104|l . Up 
to the second order, we have 



2 



B 



= t&A+B + y {^A+B^ + 5a^B - Sb^a) , (105) 



1 ~ t'SA+B, (106) 



and hence 



{e^^^e'^'^ - if {A + e'^^B) ~ {A + B) +t6AB + j5a^B, (107) 



1) {A + e''^B) ~ t5A+B[A + B) 



. 2 (-^a+b' + 5a5b - 5b5a) {A + B) + {5a + 5b) 5aB 

= -{5A6BA~5B6AB) + t'{5A + 5B)5AB, (108) 
{Qt&A^tSu [A + e'^'^B) ~ t^5A+B^ {A + B)^{). (109) 

Summing up the second-order terms with the coefhcicnt (— l)'"'/(fc + 1), we have 

^^5a^B - ^ {5a5bA - 5b5aB) -^^{5a + 5b) 5aB = ^ {5a^B + 5b^A) , (110) 
which we integrate to obtain 

^ {5a'B + 5b' A) = ^ {[A, [A, B]] + [[A, B],B]). (Ill) 
5.5 Example: Ruth's formula 

We now extend the above computation to the exponential product 

QPixA^P2xB^paxA^P4,xB^p^xA^pexB _ g$(a;) (112) 
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and seek Ruth's formula H64|) as a specific solution of the general formula. We compute the second-order 
and third-order correction terms of defined in 



$(a;) = x{A + B)+ x^R2 + x^Rs + O(x^), 

and put R2 = R3 = 0. 

The same computation as from Eq. H99|) through Eq. 11U4|) produces 



(113) 



t^A pP2tSB pPatSA iJPitSB c,P5tSA gJl^tSB 



-1)' 



{piA 



P2B + C 



.Pits A ^P2tSi 



P^A---) dt. 



(114) 



Note again that all the terms are given by commutation relations. 
For the term fc = 0, we have up to the second order of x, 

PxA + eP''^^p2B + eP'^^^'eP^'^^'psA-- ■ 

2r 2 



PlA + (^1 + tPi^A + —Pi^Sa^J P2B 

1 + t {piSa + P25b) + Y {pi'Sa^ + 2piP2SaSb + P2^Sb^) P3A -\ 

{Pl +P3+P5)A+ {p2 +P4+ Pe) B 

+t [PiP25aB + P2P3SBA + {Pl +P3)Pi^AB + (jP2 +P4)P5^S^ + {Pl + P3 + Pb) Pg^aB] 

e 
+ 2 



pIp2Sa^B + pIpsSb^A + 2piP2P3SaSbA + {pi + ps) p^Sa^B 
+ 2p2P3P4Sb6aB + {p2 + pa^psSb^A + 2 {pip2 +PiPa + P3P4) Ps^aSbA 



(115) 



+ (Pi +P3+P5) PqSa^B + 2 {p2P3 + P2P5 + PiPb) PqSbSaB . 
The zeroth-order term with respect to t appears only here and hence we have the conditions 

Pi+P3+P5 = ^ and P2+P4+P6==l- (116) 
Using Eq. (|116|l and the identity 5b A ~ —5aB, we can reduce the right-hand side of Eq. pi5|l as 



with 



A + B + t{l- 2q)6AB + Y [(1 - g - Sr)SA^B + {q - 3s)Sb'^A] 



q = P2P3+P2P5+PiP5, 

r = P1P2P3 + P1P2P5 + PlPAPb + P3PAPb , 

S = P2P3PA + P2P3P6 + P2P5P6 + P4P5P6 ■ 



(117) 



(118) 
(119) 
(120) 



For fc = 1, we first have 



oPI^Sa „P2tSBcJ>3tSAcyP4tSB„P5tSA„P6tSE 



1 ~ tSA+B + -77 [Sa^ + fe' + 2(1 - q)SA6B + 2(7^5^] , (121) 



where we already used the conditions in Eq. pi()|) . Applying Eq. 11211) to Eq. (|117|l and dropping 
higher-order terms, we note that the first-order term vanishes and have 



pPltSAcJ>2tSB aPat^AfJIitSB cJ>5tSAfJ)6tSB 



,g^...-^g^...-„g^o^-^g^.^-og^...-^g^„..-„ _ 2) [piA + ePi"-*p2-B + ePi*''^eP^"^p3A • • • ) 

<2(1 - 2q)SA+BSAB + J [Sa^ + 5b^ + 2(1 - q)5A6B + 2qSBSA] (A + B) 
^{1 - 2q) {Sa^B - Sb^A) + ^ [6a^B + Sb^A - 2(1 - q)5A^B - 2qSB^A] 



( 2 - 9 ) {Sa'B - 5b' A) 



(122) 
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The second-order term of t in the term k = 2 vanishes just as in Eq. H1U9|) . Thus we arrive at 



<P{x) =x{A + B) + -il- 2q)5AB + - 



1 



0{x^). (123) 



Putting the second-order and third-order terms to zero, we have a set of simultaneous equations of the 
parameters as 

P1+P3+P5 = 1, (124) 
P2+P4+P6 = 1, (125) 
2q = 2 {p2Pz + P2Pb + PAPb) = 1, (126) 

6r ™ 3 (pi + 2p3P4P5) = 1, (127) 

65*^3(2^2^3^4+^6) = 1- (128) 

We can confirm that Ruth's formula H64() . or 

7 2 3 2 1 

Pi = ^2^3' ^3^4' ^"^^"3' ^^"~24' ^"^^ ^^6 = 1 (129) 

is indeed a solution of the above set of simultaneous equations. With six variables for five equations, the 
solution is in fact a continuous line; Ruth's solution (|129|l is just a point on the line. By adjusting the 
last variable pe, we have the continuous solution shown in Fig.l^l (We can solve the set of equations with 
five parameters by putting p^, ~ Q, but the solution is complex.) 

5.6 Example: perturbational composition 

We finally present an interesting exercise, motivated by the "perturbational composition" [44]. Suppose 
that we apply a weak transverse field to an Ising spin. We ask what is the correction term in the exponent 
of the right-hand side of 

gf 70-xg2;CT,gf7(j^ ^ g*(a;,7) ^ ga:(o-z +7C1 (x)+0(7^)) _ {'^i^) 

Notice that we expand the exponent with respect to the perturbation parameter 7, not with respect to x 
as in the preceding sections. The first-order perturbation term Ci{x) in turn contains higher orders of x. 
We could explicitly compute the 2x2 matrices on both sides of Eq. (|130|l . expand them with respect to 7 
and compare them term by term, but the quantum analysis provides a more elegant way of computation. 
We differentiate the both sides of Eq. H130() with respect to 7: 

-^gfT'T^g^'^^gfTO-^ 

d7 



Equating the both sides, we have 
07 



de* 


a$(a:,7) _^*1 




d$" ' 


^7 


(5$ 97 ' 


X , 


gf 70-xga;cr,gf7<Tx _ 


^_gf7'Txg:!:<T,gf7a,^^ 




e (Ta;e -1- (Ta; j = 


- |e* (e-^^ + 1) a. 


) = 


" (1 + 

2 1 - e-'5* ^ 





(132) 



|(a:<5..-t-0(7))[^^'T.. (133) 



Putting 7 = 0, we have 



( \ If 1 + exp (-X^o^,) Iv^ nr n /-,o^\ 

2 1 — exp (— xOcr^) 2^-^ 



n=0 
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P, X 

Figure 5: The solution line of the Figure 6: The coefficient o f the first- 
set of simultaneous equations (|124|l - order perturbation of Eq. (|141|l . 

(inHi- 



owing to the fact (5$ — x6a^ + 0(7), where we have made the Taylor expansion 

1 + e- 



X- 



- = ^a„x" (135) 



1-e 

ri=0 

with — 1. We also note that the function H135|l is even with respect to x and hence a„ — for odd 
integers n. 

The right-hand side of Eq. (|134|l is explicitly calculated as follows. In each order, we have 

= [az,<Jx] ^2iay, (136) 
Sa, (Tx = 2i [o-^, = 40-2,, (137) 
Scr.^o-x = 'i[az,ax] = Siffy, (138) 



or in general, 



V. = i (139) 
" ' for even n. ^ ' 



We substitute Eq. (|139|l for each even-order term of the right-hand side of Eq. (|134|l and arrive at 

Ci{x) = - 2_] an(2a;)"cr3; = x-^——:^ax = (a; coth a;) cr^, . (140) 

n=0 

(In the second equality, we used the Taylor expansion (|135|l in the reverse direction.) 
In summary, we have 

gf 70-xg2;(T^gf 7(7:^ _ ^x[a^+l(xcot\ix)(T^+0(l'^)\ ^ (141) 

The coefficient x coth x behaves as shown in Fig. |B| We have x coth a; ~ 1 for small x as is expected, but 
xcotha; ~ x for large x, and hence the first-order perturbation term grows as x^ . 

6 Summary 

In the present article, we have reviewed a continual effort on generalization of the Trotter formula to 
higher-order exponential product formulas. As was emphasized in Sect. [2 the exponential product 
formula is a good and useful approximant, particularly because it conserves important symmetries of the 
system dynamics. 

We focused on two algorithms of constructing higher-order exponential product formulas. The first is 
the fractal decomposition, where we construct higher-order formulas recursively. The second is to make 
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use of the quantum analysis, where we compute higher-order correction terms directly. As interludes, 
we also have described the decomposition of symplectic integrators, the approximation of time-ordered 
exponentials, and the perturbational composition. It is our hope that the readers find the present article 
a useful and tutorial "manual" when they numerically investigate dynamical systems. For more practical 
applications of the exponential product formulas, we refer the readers to the review articles found in 
Refs. [45-48]. 

A Hybrid exponential product formula 

We mention here another kind of the exponential product formula [20-22] . Consider the Trotter approx- 
imant 

We can cancel out the second-order correction term in the form 

qxA^xB^-^x^IA.B] ^ ^xiA+B) + Oix'') ^ ^^^g-j 

If, in some problems, the commutation relation [A, B] is easily diagonalized, Eq. H143|) may be a useful 
approximant. 

A more complicated one is the fourth-order approximant [20-22] 

e^[5.[A.S]]^^ S, (I) Sa (I) eAmA,B]] ^ ^.(A+B)+0(.=)^ (144) 

where 

S-aCa;) = e^^V^e*^^ and ^'^(a;) ee e^^^e^^e^^^. (145) 
In fact, the diffusion equation is described by 

A = -^A and B^V{q) (146) 

and we have 

[B,[A,B]]^{VV{q)f. (147) 

The above type of the exponential product formula was referred to as the hybrid exponential product 
formula. We do not give its details in this article, since commutation relations are not easily diagonalized 
except for a few specific problems. 

B World-line quantum Monte Carlo method 

In the present appendix, we give a short review of the world-line quantum Monte Carlo method. The 
world- line quantum Monte Carlo method is to transform the partition function IjlOl) of a quantum system 
Ti. into the partition function of a classical system by means of the path-integral representation and 
simulate the latter system. We explain the method using the transverse Ising model (jllf) . or H = A + B 
with 

A^-J2 Jv'^i'^j and B = -T^af. (148) 

(hi) * 

The starting point is the Trotter decomposition H21|) of the partition function, namely the Suzuki- 
Trotter transformation [3], of the form: 

Z . TVe--. limTv(e-^-e-^-)"^^limE(K^^^ 

= E ({.r}|e-^-e-^-e-^-e-^-...e^^-|{.f)}). (149) 

In the second line, we have taken the trace with respect to a complete basis set by using the spin z axis 
as the quantization axis: 
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where the eigenvakie is cr^"-' = ±1. The meaning of the superscript (0) becomes self-evident just below. 
We further insert the resolution of unity in between each pair of the exponential operators in the last line 
of Eq. (|149(l . obtaining 



E E E - E ({-.'°'} 
(,.»}{.■■■){.:■>} {.:->} 



{.."'} h-1K'})(K'} 



K'}) -({..'■-«} 



e " 



{.;»'}). (151) 



In the above expression, we used the fact that the operator A is diagonal in the representation of {cr'™''} 
and hence made the complete set on the both sides of each operator e^^^ identical. In contrast, the 
operator e " has off-diagonal elements. 

Let us calculate the matrix elements in Eq. H151|) . The matrix element of the operator e~^^ is easy: 



f {rn)\\ / /3 7 (m) ("i) 1 



(152) 



This is because the operators {erf} are all diagonal in the present representation as in Eq. (|150|l . On the 
other hand, the operator e~~^ has off-diagonal elements as well as diagonal ones in the following form: 



with each matrix element given by 



(m+l) 



(m+1) 



(m) 



(m) 



-1 



V 



+^) = +l) 


^(m+l) ^ 




cosh — 
n 


sinh ^— 
n 


->> 


sinh ^— 

n 


cosh — 
n 





These matrix elements are expressed in a single equation 

/ (m) ^ct" (m+l)\ / (rn) (m+1) , . \ 

where the parameters 7„ and (5„ are defined in 



g7n+5n (2osh and e '>'"+'''» =sinh^— , 

n n 



or more explicitly defined by 

In = 



1 1 . 1 1 . 2/3r 

- - log tanh — and o„ = - log - smh . 

2 n 2 2 n 



The expressions (I152|l and (|155|l give the partition function (|151|l in the form [3] 

Z= hm y e-^^- 



(153) 



(154) 



(155) 



(156) 



(157) 



(158) 



with the resulting classical Hamiltonian [3] 



m— 2 



(159) 



(") - ^(0) 



where we dropped a constant term due to S„. Note that the periodic boundary conditions, a\ = a. 
must be required in the second term of Eq. ()159|l because the trace operation in Eq. H149|l demands it. 
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Figure 7: The three-dimensional 
classical system (|159ll mapped from 
the two-dimensional quantum sys- 
tem (PH| . 
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Figure 8: In the Trotter limit n — > cx), the Trot- 
ter axis becomes a continuum. The intra-layer 
interaction becomes an interaction between two 
continuum axes. 
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Figure 9: Spins on lattice points 
become domains on Trotter axes in 
the Trotter limit n ^ oo. 



The classical Hamiltonian (|159() is interpreted as follows (Fig.[7J. Suppose that the original quantum 
system (|148|) is defined on a square lattice. The first term of Eq. H159|) indicates that the two-dimensional 
system is replicated into n layers with the intra-layer interaction reduced by n times. The second term 
of Eq. (|159|l represents the inter-layer interactions. The coupling is —^n/fi as defined in Eq. H159|l . Thus 
the quantum system on a square lattice is mapped to an Ising model on a cubic lattice. In general, a 
(i-dimensional quantum system is mapped to a (c?-|- l)-dimensional classical system. The additional axis 
is called the Trotter direction. The physical quantities of the quantum system can be estimated by Monte 
Carlo simulation of the mapped classical system. This is the basic idea of the world-line quantum Monte 
Carlo method [3] . 

We can use this mapping in order to study the quantum annealing [40-42] . Suppose that we look for 
the ground state of the diagonal part A of the system (|148f) . Random exchange interactions { Jy} may 
produce many local minima that are only slightly above the ground state in energy but far apart from 
the ground state in the phase space. The simulated annealing, a well-known method of ground-state 
search, is often trapped in a local minima and does not reach the ground state. In quantum annealing, 
we use the transverse field F in order to induce tunneling from local minima to the ground state. We 
first apply the off-diagonal part B of Eq. H148|) strongly and turn it off gradually, hoping to end up with 
the ground state of the diagonal part A. This corresponds to a Monte Carlo simulation of the mapped 
classical system (|159|) with the intra-layer coupling 7„ being infinitesimally weak at the beginning and 
infinitely strong at the end. Each layer of the system (|159|) is first independent of each other and is 
gradually frozen into an identical configuration, which we hope is the ground state. 
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An annoying problem inherent in the algorithm of the quantum Monte Carlo method is the systematic 
error due to the finite Trotter number n. It used to be that simulations were carried out for various finite 
values of n, quantities were estimated in each simulation, and then the limit n ^ oo was taken in the 
process of the data analysis, which was called the Trotter extrapolation. A recent development of the 
quantum Monte Carlo method dramatically changed the situation. We here mention the development 
briefly; see Ref. [39] for a tutorial and exhaustive review of the topic. 

In the most recent quantum Monte Carlo algorithm, it is possible for some systems to take the Trotter 
limit before we set up the classical system for simulation. Taking the Trotter limit n ^ oo, we have a 
continuum Trotter axis (Fig. (Note again that the boundary conditions are required in the Trotter 
direction.) The interaction is described as follows (Fig. O. Instead of Ising spins on lattice points of 
a Trotter axis, we have up-spin domains and down-spin domains on the axis. Instead of intra-layer 
interactions between a pair of nearest-neighbor spins, we have parallel-spin areas and anti-parallel-spin 
areas. In Monte Carlo simulation, we update the up-spin domains and down-spin domains on the basis 
of the energy of the parallel-spin areas and anti-parallel-spin areas. 

It is thus possible in such situations to carry out a simulation in the Trotter limit n — )■ oo. Monte 
Carlo estimates of such a simulation are free of the systematic error of the order 0^ /n in Eq. H21|) . and 
hence do not need the higher-order exponential product formula for such systems. 
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